*****************************************************************
* Replication directory for                                   ***
* Prime locations                                             ***
* by Gabriel M. Ahlfeldt, Thilo N.H. Albers, Kristian Behrens ***
* Published in American Economic Review: Insights             ***
*****************************************************************
* 01/2025
* Stata
version 17.0

* This do file analyese variation in the spatial structure of cities

* Load data
	u "$temp/metro_using.dta", clear
* Merge covariates
	merge 1:1 metro_id using "$data_125cities/METRO_LEVEL_COVARIATES/metro_covariates.dta"
	drop _m
	
* Number of prime locations
	sum Number_PL
	display r(sum)
* Share relatively monocentric cities
	sum Number_PL if Number_PL == 1 & PL_emp_share > 0.25
	local share = r(N)/125 
	display `share'	

* Define simple program for easy illustration *
	capture program drop SCATTER	
	program define SCATTER	
		local ylab: variable label `1'
		reg `1' `2' , robust
			local b_all = round(_b[`2'],0.01)
			local b_all_ft : display %3.2f `b_all'
			local se_all = round(_se[`2'],0.01)
			local se_all_ft : display %3.2f `se_all'
		reg `1' `2' if  world_region =="North America", robust
			local b_small = round(_b[`2'],0.01)
			local b_small_ft : display %3.2f `b_small'
			local se_small = round(_se[`2'],0.01)
			local se_small_ft : display %3.2f `se_small'
		reg `1' `2' if   world_region !="North America", robust
			local b_large = round(_b[`2'],0.01)
			local b_large_ft : display %3.2f `b_large'
			local se_large = round(_se[`2'],0.01)
			local se_large_ft : display %3.2f `se_large'
		twoway 	(scatter `1' `2' if  world_region =="North America", mcolor(blue%25) mlcolor(blue%0) msize(small))  ///
				(lfit `1' `2' if  world_region =="North America", lcolor(blue) lpattern(shortdash)) ///
				(scatter `1' `2' if  world_region !="North America", mcolor(red%25) mlcolor(red%0) msize(small) msymbol(S)) ///
				(lfit `1' `2' if  world_region !="North America", lcolor(red) lpattern(longdash)) ///
				, legend(order(1 "North America (NA)" 3 "Rest of world (RoW)") pos(6) cols(2) region(lstyle(none))) name(`3', replace) ytitle("`ylab'") `4' `5' `6' `7' `8' ///
				note("dy/dx = `b_all_ft' (`se_all_ft')" "dy/dx | NA = `b_small_ft' (`se_small_ft')" "dy/dx | RoW = `b_large_ft' (`se_large_ft')") graphregion(color(white))
	end
* Programme definition ends *
	
* Gen log of metro pop	
	capture gen lmetropop_2015 = ln(metropop_2015)
	label var lmetropop_2015 "Ln metro population"
* Retreive maximum number of PLs for best axis labels
	sum Number_PL
	local PLnum = r(max)
* Use programme to illustrate correlations
	SCATTER Number_PL lmetropop_2015	"num_pop"  "xlabel(13[1]17)" "ylabel(1[1]`PLnum')"		 
	SCATTER PL_emp_pcshare lmetropop_2015	 "share_pop" "xlabel(13[1]17)"
	SCATTER PL_emp_pcshare Number_PL		 "share_num" "xlabel(1[1]`PLnum')"
	grc1leg num_pop share_pop share_num , graphregion(color(white)) scale(1.25) cols(3)  xsize(12) ysize(4) graphregion(margin(0 0 0 0)) plotregion(color(white))
* Export Figure 2, panel (b) ***************************************************	
	capture mkdir "$figures/GlobalCities"
	graph export "$figures/GlobalCities/FIG_2b_GlobalCities-correlations.png", replace  height(600) width(1800) 
* drop log metro pop variable
	capture drop lmetropop_2015	

* Compute Number and importance of prime locations by world region	
	preserve
	* Aggregate to world region
		replace world_region = "Africa & Middle East" if world_region == "Africa" | world_region == "Middle East"
		collapse (count) CityN=Number_PL (mean) Number_PL PL_area_share PL_emp_share  , by(world_region)
		foreach var of varlist PL_area_share PL_emp_share  {
			replace `var'=`var'*100
			format `var' %9.2f
			tostring `var', replace force format(%9.2f)
			replace `var'=`var'+"%"
		}			
		format Number_PL %9.2f
		tostring Number_PL, replace force format(%9.2f)
		label var world_region "Region"
		label var CityN "Cities"
		label var Number "PLs per city"
		label var PL_area_share "PL share at area"
		label var PL_emp_share "PL share at emp"
		capture mkdir "$tables_App/GlobalCities"
* Export Table B.4.1 ***********************************************************
		texsave * using "$tables_App/GlobalCities/TAB_B4_1_GlobalSummaryStats.tex", title("Number and importance of prime locations by world region") ///
		size("footnotesize")   width(15cm) align(lccccc)  varlabels replace  frag  ///
		footnote("Agglomerated and dispersed cities are defined by on whether the share of prime locations at total employment is larger or smaller than the median share.")	
	restore


 * Variable generation for further analysis
	gen lPL_emp_share = ln(PL_emp_share+0.01)
	gen lNumber_PL = ln( Number_PL)
	gen subway = OPEN_YEAR <=2015
	gen SubwayYears = 2015-OPEN_YEAR 
	replace SubwayYears = 0 if SubwayYears == .
	foreach var of varlist  subway metropop_2015 sWATER sSTEEPSLOPE{ 
		gen l`var' = ln(`var')
		sum l`var'
		replace l`var' = l`var' - r(mean)
		sum `var'
		replace `var' = `var'-r(mean)
		}
	encode world_region, gen(world_region_id)
	encode country, gen(country_id)
	tab world_region, gen(WRD)
	gen WR_AfricaME=1 if world_region=="Africa" | world_region=="Middle East"
	gen WR_Asia=1 if world_region=="Asia"  
	gen WR_Australia=1 if world_region=="Australia"  
	gen WR_Europe=1 if world_region=="Europe"  
	gen WR_SouthAmerica=1 if world_region=="South America"  
	foreach x in  WR_AfricaME WR_Europe WR_Asia WR_Australia WR_SouthAmerica{
		replace `x'= 0 if `x'==.
		}

* Is the negative relationship between PL share and number of PLs robust?	
	label var Number_PL "\# prime locations"
	label var lPL_emp_share "Ln employment share"
	label var lNumber_PL "Ln \# prime locations"
	label var lmetropop_2015 "Ln metro population"
	label var subway "Subway"
	label var SubwayYears "Years since subway opening"
	label var  stations_2010_TA_edit "\# subway stations"
	label var sWATER "Share water (\%)"
	label var  sSTEEPSLOP "Share steep slope (\%)"
	
	eststo: reg PL_emp_pcshare Number_PL , robust	 
		estadd local FE = "-"	
	eststo: reg PL_emp_pcshare Number_PL lmetropop_2015 sWATER sSTEEPSLOPE  subway  SubwayYears stations_2010_TA_edit , robust 
		estadd local FE = "-"		
	eststo: reg PL_emp_pcshare Number_PL lmetropop_2015 sWATER sSTEEPSLOPE  subway  SubwayYears stations_2010_TA_edit, robust  abs(country_id)
		estadd local FE = "Yes"	
* Export regression results - results discussed in text at the end of Section 3 
	capture mkdir "$tables_App/GlobalCities"	
	esttab  using "$tables_App/GlobalCities/FIG_B4_Subadditivity_RobustnessTable.tex", replace b(3) se(2) label compress  wide r2(3) stats(FE N r2 , fmt(%18.3g )   ///
	labels("`Region effects'"  `"Observations"' `"\(R^{2}\)"' ))  ///   stats(Country_effects IV N r2) 0.* 
	title("Determinants of prime locations employment share") modelwidth(6) nogap star(+ 0.15 * 0.1 ** 0.05 *** 0.01) ///
	addnote( "Unit of observation grid cell. Undevelopable grid cells excluded. PL is a dummy taking the value of one if a grid cell belongs to at prime location, and zero otherwise." )	
			eststo clear	
	
* Generate Table 2: Determinants of the geography of prime locations
	
	* Panel B
	eststo: reg PL_emp_pcshare WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica , robust	
		estadd local BaseControls = "-"		
		estadd local SubwayControls = "-"			
	eststo: reg PL_emp_pcshare  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica lmetropop_2015 sWATER sSTEEPSLOPE  , robust
		estadd local BaseControls = "Yes"		
		estadd local SubwayControls = "-"			
	eststo: reg PL_emp_pcshare  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica lmetropop_2015 sWATER sSTEEPSLOPE   subway SubwayYears , robust	 
		estadd local BaseControls = "Yes"		
		estadd local SubwayControls = "Yes"		
	eststo: reg PL_emp_pcshare  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica lmetropop_2015 sWATER sSTEEPSLOPE   subway SubwayYears stations_2010_TA_edit, robust	 
		estadd local BaseControls = "Yes"		
		estadd local SubwayControls = "Yes"	
* Export Table 2, panel B ******************************************************
	capture mkdir "$tables/GlobalCities"	
	esttab  using "$tables/GlobalCities/TAB_2b_PanelB.tex", replace b(3) se(3) label compress  wide r2(3) stats(BaseControls SubwayControls  N r2 , fmt(%18.3g ) ///
	labels("`Basic controls'" "`Subway controls'"  `"Observations"' `"\(R^{2}\)"' ))  ///   stats(Country_effects IV N r2) 0.* 
	title("Determinants of prime locations employment share") modelwidth(6) nogap nostar addnote( "Unit of observation grid cell. Undevelopable grid cells excluded. PL is a dummy taking the value of one if a grid cell belongs to at prime location, and zero otherwise." )	// star(+ 0.15 * 0.1 ** 0.05 *** 0.01)  
			eststo clear
		
		* Panel A
	tab Number_PL if PL_emp_share > 0.25	
	eststo: reg Number_PL  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica, robust	
		estadd local BaseControls = "-"		
		estadd local SubwayControls = "-"
	eststo: reg Number_PL  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica  lmetropop_2015 sWATER sSTEEPSLOPE  , robust	
		estadd local BaseControls = "-"		
		estadd local SubwayControls = "-"
		eststo: reg Number_PL  WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica  lmetropop_2015 sWATER sSTEEPSLOPE   subway SubwayYears, robust	
		estadd local BaseControls = "-"		
		estadd local SubwayControls = "-"
		eststo: reg  Number_PL WR_AfricaME  WR_Asia WR_Australia WR_Europe WR_SouthAmerica lmetropop_2015 sWATER sSTEEPSLOPE   subway SubwayYears stations_2010_TA_edit, robust	
		estadd local BaseControls = "-"		
		estadd local SubwayControls = "-"
* Export Table 2, panel A ******************************************************			
	esttab  using "$tables/GlobalCities/TAB_2a_PanelA.tex", replace b(3) se(3) label compress  wide r2(3) stats(BaseControls SubwayControls  N r2 , fmt(%18.3g )   labels("`Basic controls'" "`Subway controls'"  `"Observations"' `"\(R^{2}\)"' ))  ///   stats(Country_effects IV N r2) 0.* 
			title("Determinants of prime locations employment share") modelwidth(6) nogap nostar addnote( "Unit of observation grid cell. Undevelopable grid cells excluded. PL is a dummy taking the value of one if a grid cell belongs to at prime location, and zero otherwise." )	// star(+ 0.15 * 0.1 ** 0.05 *** 0.01) 
			eststo clear		
	
* Script ends	 
